% Code for Figure 3
clc
clear all
% close all

A = [0.5 -0.3;
     -0.8 0.1];
A5 = [-0.11 0.03;
      0.17 -0.11];
D = [-0.2;0.1];

nx = size(A,1);
nd = size(D,2);

H0 = 1*eye(nx); 
Hw = 0.1*eye(nd);

Abar = [A zeros(nx,nx) zeros(nx,nx) zeros(nx,nx) zeros(nx,nx) A5;
        eye(nx) zeros(nx,nx) zeros(nx,nx) zeros(nx,nx) zeros(nx,nx) zeros(nx,nx);
        zeros(nx,nx) eye(nx) zeros(nx,nx) zeros(nx,nx) zeros(nx,nx) zeros(nx,nx);
        zeros(nx,nx) zeros(nx,nx) eye(nx) zeros(nx,nx) zeros(nx,nx) zeros(nx,nx);
        zeros(nx,nx) zeros(nx,nx) zeros(nx,nx) eye(nx) zeros(nx,nx) zeros(nx,nx);
        zeros(nx,nx) zeros(nx,nx) zeros(nx,nx) zeros(nx,nx) eye(nx) zeros(nx,nx)];
Dbar = [D;zeros(5*nx,nd)];
Hbar0 = zeros(6*nx,6*nx);
Hbar0(1:nx,1:nx) = H0;
Hbark = Hbar0;

Nk = 15;
Xk = zeros(nx,1)+eye(nx)*unitbox(nx,1);
vk(:,1) = volume(Xk);

HHk = [zeros(1,nx);eye(nx)];
ck = [1;0;0];
XXk = ck+HHk*unitbox(nx,1);        
    
figure(1)
hold on
Options.color='g';
Options.alpha=0.2;
plot(XXk,Options)
    
for k=2:Nk
    Hbark= [Abar*Hbark Dbar*Hw];
    Hbark = reduceG(Hbark,10);
    
    s = size(Hbark,2);
    Xk = zeros(nx,1)+Hbark(1:nx,:)*unitbox(s,1);
    vk(:,k) = volume(Xk);
    
    HHk = [zeros(1,s);Hbark(1:nx,:)];
    ck = [k;0;0];
    XXk = ck+HHk*unitbox(s,1);        
    
    figure(1)
    hold on
    Options.color='g';
    Options.alpha=0.2;
    plot(XXk,Options)
end

figure(2)
hold on
plot(1:Nk,vk,'go')

 